A rapidly converging domain decomposition method for the 

Helmholtz equation 



Christiaan C. Stolk 

University of Amsterdam, Korteweg-de Vries Institute for Mathematics, P.O.Box 94248, 1090 GE Amsterdam, 

The Netherlands 



Abstract 

A new domain decomposition method is introduced for the heterogeneous 2-D and 3-D Helmholtz 
equations. Transmission conditions based on the perfectly matched layer (PML) are derived that 
avoid artificial reflections and match incoming and outgoing waves at the subdomain interfaces. 
We focus on a subdivision of the rectangular domain into many thin subdomains along one of the 
axes, in combination with a certain ordering for solving the subdomain problems and a GMRES 
outer iteration. This leads to very small iteration numbers in examples, independent of problem 
size and number of subdomains. When combined with multifrontal methods, the solver has near- 
linear cost. It is to our knowledge only the second method with this property next to the moving 
PML sweeping method. 
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1. Introduction 

In this paper we introduce a new domain decomposition method for the solution of the Helmholtz 
equation in two and three dimensions. To be specific we consider in 2-D 

- d 2 xx u(x, y) - d* y u(x, y) - k(x, y) 2 u(x, y) = f(x, y), (1) 

where k{x,y) — , with c{x,y) the wave speed. The computational domain is assumed to be 
a rectangle that is truncated using the perfectly matched layer pp. We focus on solving the large 
linear systems resulting from discretization with standard 5 or 7 point finite differences. 

We have two main findings. First, we have constructed new transmission conditions. These are 
designed to ensure that 

(i) the boundary conditions at the subdomain interfaces are non-reflecting; 

(ii) if Qj—i and Qj are neighboring subdomains then the outgoing wave field from 

(2) 

Clj-i equals the incoming wave field in ftj at the joint boundary and vice versa. 

This is achieved in a simple and accurate way using PML boundary layers added to the subdomains 
and single layer potentials. See [2] for a related approach in the finite element discretization of the 
time harmonic Maxwell equations. 

Our most remarkable finding concerns the situation where the domain is split into many thin 
layers along one of the axes, say J subdomains numbered from 1 to J. Following [3] we will also 
call these quasi 2-D subdomains. Generally, an increase in the number of subdomains leads to 
an increase in the number of iterations required for convergence. Here we propose and study a 
method where the number of iterations is essentially independent of the number of subdomains. 
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A necessary condition for this is that information can travel over the entire domain (or at least 
an O(l) part thereof) in one iteration. To achieve this we use a multiplicative method, where 
the subdomains are first solved consecutively from j = 1 to j = J, each time using information 
from the solution of the neighboring, previously solved subdomain, and then in the same way from 
j = J downto j = 1 using the residual as right hand side. In this way information can travel all 
over the domain with only two solves per subdomain. The procedure is used as a preconditioner 
for GMRES. 

We studied numerically the convergence of the method for different choices of the grid distance 
h and the frequency ui, keeping uh constant, and different numbers of subdomains. In our examples 
the method converged rapidly, with generally less than 10 iterations needed for reduction of the 
residual by 10~ 6 . Moreover, the required number of iterations was essentially independent of the 
size of the domain and the number of subdomains. 

This is attractive in combination with the use of multifrontal methods for the subdomain 
solves. Indeed, as was argued by Engquist and Ying [3J, in 3-D the set of quasi 2-D subproblcms 
can be solved by multifrontal methods in 0(N log N) time, with 0(N 4 / 3 ) cost for the factorization, 
versus 0(N 3 ^ 2 ) and 0(N 2 ) when the multifrontal method is applied directly to the 3-D system. 
The method therefore behaves near-linearljQ It is the second such method we are aware of, in 
addition to the moving PML sweeping method, for which such observations were made in [3J. 

Two things have to be kept in mind. Of course the method is only near-linear provided that 
solutions are required for a sufficiently large number of right hand sides to recoup the cost of the 
factorization. Secondly we estimate, based on our examples, that the thickness of the PML layers 
needs to increase with increasing N. A required growth of O(logiV) is consistent with our data. 
This would lead to an additional factor 0(log N) for the cost of the solves and 0((\ogN) 2 ) for the 
cost of the preparation. 

1.1. The method and its context 

Next we discuss in more detail the ideas behind the method and some of the relevant literature. 

To motivate our approach we recall the 1-D problem with k = constant, see the review in [4] or 
[HE! [7]. Let ]0, L[ be the domain. The differential equation and the Robin boundary conditions 
read 

—dl x u(x) - k 2 u(x) = f(x) for < x < L, 

d x u + iku = at x = 

—d x u + iku = at x = L. 

The Robin boundary conditions are exact non-reflecting boundary conditions and ensure that 
there are no incoming waves at the boundaries. We assume the domain is divided in J subdomains 
]bj-i, bj [ with 

= b a < bi < . . . < bj = L. 

The original problem is then equivalent to J subdomain problems with continuity conditions at 
the interfaces as follows 

-d 2 xx u {j) - k 2 u {3) = f& for x e]6j_i, 6j [ 

d x u (j) + iku {3) = d x u {j - 1] + iku {3 - 1] at x = bj-i 

-d x u {3) + ikuW = - d x u {3+1) + iku {j+1) &tx = bj 

(by convention = = u^ J+1 ^). These continuity conditions satisfy the property To obtain 
an iterative solution method, the right hand side of the continuity conditions is taken from the 
previous iteration, i.e. a sequence v rL is constructed, where n is the iteration number and j the 
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subdomain index according to 

-d 2 xx v^ - k 2 v%) = f U) fo,. x € ] bj _ 1)bj [ ( 3 ) 

d x v$ + ikv& = d x v^ ] + ikv ( Jlp , at x = (4) 

^ + ifct#> = - a^-i + ^wi J -i x) at * = &j ■ ( 5 ) 

This method is optimal in the sense that it converges in a finite number, namely </, of iterations. 
Indeed, recall that the solution for the problem —d xx u(x) — k 2 u(x) = f{x) with Robin boundary 
conditions d x ii(0) + iku(0) = hi, —d x ii(L) + iku(L) = hi is given by 

It follows by induction, starting from Vq ^ = 0, that viP satisfies 



^ } (x) = ^\ a e^ x - s) f(s) ds +2kj x e- mx - s) f(s) ds 

where A — 6 ma x(oj-n) an d B — 6 m in(,/j+n-i)- After J steps, ^4 = and B = L for all j G 
{1 I\- 

This work answers two questions about the iterative method (|3j)- ([sj . The first question concerns 
the generalization of the transmission conditions to two and three dimensions. The Robin boundary 
conditions then no longer satisfies the properties ^ (see the argument around (24) below). Several 



approaches have been proposed in the literature. First, the Robin boundary conditions can still 
be used as transmission conditions [8j. Several authors have also considered optimized Robin 
transmission conditions IIP] . A second possible approach follows from work about absorbing 



boundary conditions, e.g. [TT]. Pade approximations for A in (27) (see below) can be used to obtain 
numerical absorbing boundary and transmission conditions [121 113) . In this paper we use PML 
boundary layers [T] to achieve ^ . See [HI [2] for earlier work using domain decomposition with 
PML's. 

The second question concerns the case of large J. In one iteration of ([3|-([5]) information 
from one subdomain can only travel to its neighbors. The method therefore requires at least 
O(J) iterations to converge, hence 0(J 2 ) subdomain solves. On the other hand, by using the 
multiplicative approach outlined below, solving the subdomains consecutively, first j = 1, 2, . . . , J, 
and then j = J, J — 1, .... 1, information can travel over the full domain in just 2 solves per 
subdomain. Here we will follow this multiplicative approach. 

As mentioned out, the case of a large number of thin layers, say k grid points thick, is of interest 
when the method is used in combination with multifrontal methods for the subdomain solves. The 
computational cost of such a setup was analyzed by Engquist and Ying [3], using results of [15] , 
Consider a cube with n x n x n gridpoints, hence N = n 3 . The cost of a LU decomposition of a 
subdomain of the form n x n x k is 0(fc 3 n 3 ), while the cost of a backsubstition is 0(k 2 n 2 logn). 
Assuming k = O(l), the total cost of the factorizations is 0(7V 4 / 3 ), while the total cost per iteration 
is 0(N log N). If the number of iterations depends weakly on problem size, as we see in examples, 
then this method scales well. In the presence of PML layers that have a thickness of w pm \ grid 
points, the value k = 2u> pm i is the optimal value for the thickness of the PML layers, i.e. minimizes 
the cost of applying one set of subdomain solves. The details of our method will be explained in 
section [2] 

The papers [TBI H7] provide a review of solution methods for the Hclmholtz equation. Some 
recent other work is given in [T51 [TH] . 

1.2. Results 

Our first main result is a theoretical result, concerning the constant coefficient problem on 
a strip. Assuming that the PML layers perfectly reproduce the behavior of the solution on the 
unbounded domain, the methods solves this problem in one iteration, i.e. in one upward and one 
downward sequence of solves. We observe that the upward and downward sequence of solves can in 
fact be performed simultaneously, if the point where the the sequences cross is handled carefully. 

The second main result is the good convergence behavior in numerical examples that was 
already mentioned in the first part of this introduction. 
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1.3. Contents 

The paper is organized as follows. The next section explains in detail our method. Section [3] 
contains some theoretical results. Then in section [4] the numerical examples are discussed. We end 
the paper with a short discussion. 

2. The method 

2.1. Continuous formulation 

In this section we formulate our method in 2-D. The domain is assumed to be a set of the form 
O =]0, L[x]0, 1[. It is straightforward to generalize this to rectangular domains of different size, 
and to 3-D rectangular domains. 

The Hclmholtz operator will be referred to as A, given away from the PML boundary layers by 

A=~d 2 xx -d 2 yy -k(x,y) 2 . 
The operator in a PML layer at a boundary, say x = constant, is obtained by replacing 

d_ 1 d 

dx i + i°^Mdx 

UJ 

where a x = in the interior of the domain, and positive inside the PML layers [30] . 

The domain is divided into J subdomains along the x-axis. The interface locations will be 
denoted by x — bj , where 

= b < ■ ■ ■ < bj = L 

The "core" subdomains, without additional PML layers, will be denoted by D^> =]bj—i, bj[x]0, 1[. 
With PML layers added the notation fiO) will be used. The latter sets are obtained by padding 
the D<J) with PML layers of size Lpmi at the internal boundaries, i.e. 

Q U) =}bj-i - Lpml(1 - S jA ),bj + L PML (1 - ^,j)[x]0,l[ 

On the domains fl^\ functions few (x, y) are defined that agree with k on , and are independent 
of x and equal to k at the boundary of the core subdomain inside the added PML layers, i.e. 

{k(x,y) for bj_i < x <bj 
fc(6j-i, y) for x < bj-! (if j > 1) 
k(bj,y) for x > b 3 (if j < J). 

On the domains Cl^' operators are defined as Helmholtz operators with PML modifications, 
similar as A was defined on f2. 

Next we consider the approximation by domain decomposition of a solution u to the 2-D 
Hclmholtz equation Au = f. The function / is assumed to be integrable, which allows the definition 
of on fiW by 

(i) = f f(x) iixG DO') 
\ otherwise 

A first set of subdomain solutions ti«' is obtained by solving the equations 

A U) V U) = f U) _ 2S(x - bj^dvvV-Vibj-u ■), (7) 

consecutively for j = 1, 2, . . . , J. Here by convention = 0. A function v on Q is then defined 

by 

v(x,y) — v^'(x,y) with j s.t. bj-\ < x < bj. (8) 

The second term in the right hand side of Q requires some explanation. While this is mostly 
done in the next section, a short intuitive explanation goes as follows. The term w^ -1 ) (b .) 
exclusively contains forward going waves because of the presence of a PML non- reflecting layer 
immediately to its right in fiw). The term — 25(x — 6j_i)9 x ?;0' _1 )(6 J ._ 1) .) i s meant to cause the 
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same forward going wave field in the field v^K The form of this term can be explained by the 
properties of the single layer potential. The solution to 

Au = h(y)S(x — bj-i) 

has the property that 

limd x u(&,-_i +e,y)- d x u(bj-i - e,y) = -h(y). 

e— >0 

Assuming the medium is independent of x, the source h(y)5{x — bj_i) generates waves propa- 
gating both forwardly and backwardly in a symmetric fashion. The factor —2 is introduced so that 
the forward propagating part equals y). The backward propagating part is absorbed 

in the neighboring PML layer. Note that in this way, all the subdomain sources /( fe ) with k < j 
can contribute to the field v^K 

The downward sequence of subdomain solves takes as right hand side the restrictions to a 
subdomain of the residual 

g = f - Av 

However, v is undefined and generally discontinuous at the boundaries x — bj, j = 1, . . . , J — 1. 
While g is still well defined, it only exists as a generalized function (distribution), with most 
singular term of the form S'(x — bj)h(y). 

The problem with this is not that g is unsuitable as a right hand side. Solutions to Hclmholtz 
equations with distributional right hand sides in general exist. And, as a Helmholtz equation 
is formally an elliptic equation, the solutions are smooth away from the singular support of the 
right hand side. However, the restriction of g to the subdomains D^) is not well defined. Indeed, 
such a restriction is obtained by multiplying g by the indicator function Ijju) of and this 

multiplication is in general not well defined, because of the overlapping singular supports. 

Therefore we introduce a second set of domain boundaries 

= b < h < . . . < bj = L. 

with bj 7^ b k for all < j < J and < k < J. Similarly as above we defined sets and 
by D^) =]bj- u 6j[x]0,l[, and 0« =]6j_i - L PML (1 - S jA ), bj_ + L PMh (l - 5 jt j)[x]0, l[,_and we 
let be the Helmholtz operator with PML modification on fl^K The function on can 
now be defined by 

g {3 \x,y) = l£,u)g{x,y). 

Next a series of functions on fiW is determined for j = J, J — 1, . . . , 1 (computed in this 
order) from the equations 

iCO^O) = g U) + 2S(x - U j ))d x w^ j+1) (b^ ■). ( 9 ) 
and a function w is defined by 

w(x, y) — w^\x, y) 

where j is such that bj-\ < x < bj. The function w is in general undefined for x = bj, 1 < j < J — 1, 
where it is discontinuous, but this is not a problem (we don't go into detail regarding the regularity 
of the solutions in this work.) 

The approximate solution to the Helmholtz equation is given by v + w. We define P to be the 
map f i y v -\- vu . The map P can be used as a left- or right preconditioner in an iterative solution 
method like GMRES. 

2.2. Discrete formulation 

Discretization is done using finite differences. We focus on a relatively simple scheme, using 
the standard second order approximation to the Laplacian. Because the emphasis in this work is 
on the convergence of the iterative method for the discrete system, and on a proof of principle, 
questions related to the use of higher order discretizations and the use of different schemes such 
as finite elements are relegated to later work. 
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The grid distance is assumed to be equal in x and y directions and is denoted by h. The grid 
size is denoted by n x x n y . In 2-D, the finite difference approximation to Au is given by 

(Au)ij = (-v,i-i,j + 2u i:j - Uj + i,j) + (-Uij-i + 2mj - Uij+i) - kljUij 
In 3-D we use 



(Au)ij^k — To {~ u i-l,j,k + 2uij,k — Ui+l,j,k) + To (— «i,j-l,fe + 2Uij,k — Ui,j+l,k) 



1 

In the PML layers we use the approximation 

a x d a {a x d x u[Xi,yj)) = a x (x t ) 

where a x {x) = — .l x(x) ■ The subdomain boundaries are assumed to be at half grid points 
bj = xp. + i/2- The discrete equivalent to the interval ]bj—i,bj[ is therefore the set of points 

The use of two sets of subdomains, with two sets of LDL 1 factorizations of the A^> is not very 
attractive. Fortunately it is not needed. After the first set of discrete subdomain boundaries f3j is 
chosen, the second set is defined by 

A>=A> 
Pj = Pj 

&=/3,- + l for j = l,...,J-l. 

The domain for the operators A^) is given by the grid 

{xp j _ x j r \, . . . ,Xp} x {yi, . . . , y Uy } extended with PML layers of thickness u> pm i 
on the internal boundaries. 

Finally, we need to specify the derivative d x and the distribution S(x — bj) on the right hand 
side of ([7]) and (|9| We approximate derivative on a half-grid point by 

d x u(x J+1/2 ) w \{u j+1 - Uj). 



The S function is approximated by 



S(xi - x j+ i /2 ) » I 



i if + l/2)| = 1/2 

otherwise 



We generally aim that all subdomains have approximately the same size in number of gridpoints. 
Since this size is given by (n x + ( J — l)(2w pm i + l))n y , we choose the /3j such that 

Pj~Wpml+J J (11) 

2.3. Algorithm 

In the previous sections the operators A and P where specified. Our plan is to use GMRES for 
one of the following two equations, the right preconditioned system 

APv =f, u = Pv (12) 

or the left preconditioned system 

PAu = Pf (13) 



G 



Algorithm 1: Preparation 
given J, n x determine subdomain boundaries from (11) 
create matrix of operators Aw on subdomains given in ( 10 1 
perform LDL 1 decomposition 

Algorithm 2: Transform to have data only at boundaries 



solve P and u from (14) 
output <f) — f — Au 

Algorithm 3: Apply AP to an input / 
for j = 2, J, solve 
compute the residual g = / — Av 
for j = J - 1, . . . , 1, solve ^ 
compute the residual h — g — Aw 
output / — h 

Algorithm 4: Apply P to an input tp and add u 
for j = 2, J, solve Q with / = V 
compute the residual g = ip — Av 
for j = J - 1, . . . , 1, solve ([9]) 
output U + W + 

Algorithm 5: Solve Au = f 
use algorithm 2 to compute <p 

apply GMRES with AP given by algorithm 3 to solve (15) 
use algorithm 4 to compute solution u from if> 



Table 1: List of algorithms. Algorithms 1 and 5 form the top-level part of the program. 



These appear to be systems of size n x n y x n x n y , but as is common in domain decomposition 
methods, a modification of the problem to one involving only degrees of freedom near the boundary 



is possible at least for ( 13 ) 



Indeed, the GMRES iteration of the right-preconditioned problem can straightforwardly be 
restricted to the 2( J — 1) layers of grid points at Xf-j, k = f3j + 1 and k — /3y + 2, for j = 1, .... J — 1. 
This is based on two observations. The first is that for any /, the residual f — APf is only non-zero 
at grid points Xk.i with k = j3j + 1 or k = j3j + 2, for some j, 1 < j < J — 1. The second is that the 
right hand side / can easily be replaced by a right hand side <fi with the same property. Namely, 
let be restriction of / to the X/. i with k £ {Pj-i + 1, . . . , f3j}, and u^" 1 be the solution to 

^O'JfiCi) = f(J) (14) 

and set itk,i = ui] if k € {Pj-x + 1, . . . , Pj). Then as new right hand side the residual can be used 

<f> = f- Au. 

Then, after solving ip from 

Aip = 4> (15) 

using the right-preconditioned equation in the reduced space, the solution of the original problem 
is obtained by taking 

u = ip + u. 

This concludes our outline of the method. In Table [l] the main steps that can be used in a 
computer implementation are outlined. 



3. Theoretical results 

3.1. Multiplicative domain decomposition with upward and downward sweeps in 1-D 

Here we study our approach of using upward and downward sweeps of subdomain solves for 
the 1-D problem. We establish that the constant coefficient 1-D problem is solved in one step with 
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this method. A similar result holds when the upward and downward sequences of solves are done 
concurrently. 

For the upward sweep, consider v^' defined by 

-9 2 xx v U) - k 2 v^ = /W) for x €]bj_x, bj[ (16) 

d x v U) + ikv U) = d^-^ + ikv {j - x \ at x = bj-i (17) 

-d x v ij) + ikv ij) = a,tx = b j . (18) 

Then by induction it follows that 

v U) (x) = ±- f X e ik ^f(s) ds+^ [' e-^-^f(s) ds (19) 



for bj-i < x <bj. Indeed, if (191 is satisfied with j replaced by j — 1, it follows that 

d x v( j - 1) (b j - 1 )+ ifcvk'- 1 ) (&,•_!) = / ' 1 e lk ^- 1 -^f(s)ds, (20) 



which together with (|6| implies (19). 
Next let w^> satisfy 

-Cy j) - A: 2 w (i) = f U) for x G]6i-i, bj[ (21) 

d x w {i) + ifcw (i) = ^v 0-1 ) + ikv^"^, at x = (22) 

-d x w ij) + ikw [j) = - d x w { ' J+1) + ikw u+1) at x = bj, (23) 

obtained by setting = v^ J ' and solving for j = J — 1, J — 2 . . . , 1 in that order. Then by 
induction 

-<9 x w (j+1) (6 J ) + *fcw°' +1) (^) = / e- lfc(b ^ s) /(s)ds 

and for £>j_i < x < bj 

w(j) {x) = h[ elk(x ~ s) f^ ds+ ik£ e ~ ik(x ~ s) f( s ) ds 

i.e. the solution of the full problem. 

Next we discuss the case of an upward and a downward sweep that proceed concurrently. More 
precisely described it is the same method as in the introduction, but at time level n only the 
functions and Vn +1 are updated. Using the solution formula ([©J) again it can be verified 
that vj(x), given by Vj for bj-i < x < bj, is the solution of the original problem. 

3.2. PML based transmission on the strip 

Here we consider the problem with k = constant on the strip ]0, L[x]0, 1[, with Dirichlet 
boundary conditions at y = and y — 1 and PML boundary layers at x = and x = L. In 
this section we assume that a PML boundary layer behaves like a perfect non-reflecting boundary 
condition. 

The behavior of a perfect non- reflecting boundary is most easily described in the Fourier domain. 
After a Fourier transform u = J^i sin(27rZ)u;(x), 1 = 1,2,..., and writing ui(x) = u(x, r]), rj = 2wl, 
the Helmholtz equation becomes a family of ODE's that reads 

-% x u + r, 2 u-k 2 u = f(x,r ) ) (24) 

We assume that k ^ 2wl for all integers I > 0. The non-reflecting boundary condition becomes 

d x u + Xu = hi at x = (25) 

—d x u + Xu = hi at x = L, (26) 
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where A is given by 



A = 



J i\Jk 2 — r\ 2 if \rj\ < k 
-\/i] 2 — k 2 if \rj\ > k, 



(27) 



and hi and hi are for homogeneous non-reflecting boundary conditions and non-zero if incoming 
waves are to be modeled. (In the spatial domain, after inverse Fourier transform in y, the factor 
A would become a pseudodifferential operator that is non-local, explaining why in two and three 
dimension we can not obtain the properties (i) and (ii) of the introduction using Robin boundary 
conditions.) 

We have the following result: 

Theorem 1. In the situation just described, the map P satisfies APf = f. 

Proof. The solution formula for ( 24p6 ) is given by 



I 

2A 



-1 

2A 



Xx p —A(x—L) 

i-**-)f(8, V ) ds + —to + 2A to 



First we consider the fields in other words the forward sweep. Using induction, it easy to 
show that 



v^\x, rj) — 



-1 
2A 



e A(x_a) /(».»?) ds + 



-1 
2A 



,-\(x-s) 



f(s,rj) ds. 



(28) 



Indeed, assuming this is true with j — 1 substituted for j it follows that 

-I 



? x{x - s) f(s,v)ds 



The solution formula applied to right hand side f^ (x, rj) — 26 (x — bj-i)d x v^ 1 - 1 tj) then gives 
& 

Next we consider the backward sweep. The vfi' are solutions to Helmholtz equations with as 
right hand side the residual / — Av derived from the v^K From (28) it follows that 



(d x + X)v(b~j) 



= A(x-s) 



f(s)ds. 



(29) 



It follows that v + wU' satisfies for bj—i < x < bj the equations 

-d 2 xx (v + u;Cfl) + [rf - k 2 ){v + w&>) = f 
while at the boundaries of the interval 



Uh - _ 



,X(x- 



>f(s)ds 



(-d x + \){v + w^(bj - 0)) = (-d x + X)(v + w^+V) 



at x = bj—i 
at x = bj 



(30) 
(31) 



Here w^\bj — 0) denotes the limit lim^g w^\x). The first of these two equations follows easily 



Ef6 

from (29), while the second follows from the transmission condition. Then by induction (31) can 
also be written as 

(-d x + \)(v + w U) (b J -0)) = - f M 
It follows that 



f(s)ds 



at x 



v(x) + (x) 



-1 

2A 



e Mx - s) f(s,v) ds 



2A 



e -**-)f(8,T))d8 



for bj-i < x < bj, which completes the proof. 



□ 
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N x X Ny 


h (m) 




Number of x-subdomains 


3 


10 


30 


100 


300 


600 x 212 


16 


12.5 


4 


5 


6 






1175 x 400 


8 


25 


5 


6 


7 






2325 x 775 


4 


50 


6 


6 


7 


9 




4625 x 1525 


2 


100 


6 


6 


7 


8 




9225 x 3025 


1 


200 




7 


8 


9 


13 (8) (*) 



(*) 13 was obtained for w pm \ = 5, 8 for w pm i = 6. 



Table 2: Convergence results for example 1. Displayed is the number of iterations for reduction of the residual by 
10~ 6 as a function of the size of the domain and the number of subdomains. 

4. Numerical results 

In this section we present examples in 2-D and in 3-D with constant and variable k. We'll 
focus on the convergence of the method, measured by the number of iterations for reduction of 
the residual by a factor 10~ 6 . In our 2-D example we will vary the size of the domain and the 
number of subdomains, keeping hio constant. We will see that the number of iterations required 
is essentially independent of those parameters. 

In 3-D we take subdomains of constant thickness of 10 grid points, excluding the PML layers. 
The number of subdomains is therefore dictated by the size of the domain, and we study the 
convergence as a function of domain size. Again the number of iterations is approximately constant. 
We also study the influence on the parameter w pm \ for constant coefficient media. In our 3-D 
examples a value of u> pm i = 4 generally produced a good convergence. Nevertheless the parameter 
w pm i has some influence and some insight in this is obtained from the third example. The 3-D 
examples were done with domain sizes up to (400) 3 . 

Because of the size of the problems, the implementation was done under MPI. For the solution 
of the linear systems on the subdomain the parallel sparse multifrontal solver MUMPS 21j was 
used. In the version used to generate the 2-D examples the sequential sparse multifrontal solver 
UMFPACK [22] was used. The examples were run on the LISA linux cluster of the Stichting 
Academisch Rekencentrum Amsterdam (SARA). 

4-1. Example 1: Marmousi 

Our first example is the Marmousi model, a synthetic model from reflection seismology. In this 
model the velocity c(x,y) varies between 1500 and 5500 ins -1 . The model and a solution to the 
Hclmholtz equation are given in Figure [I] 

Our first set of computations shows the number of iterations required for convergence as a 
function of grid size h and the number of subdomains J. It is summarized in Table [2] The grid 
size varies between h = 1 and h = 16 m, and the number of subdomains between 3 and 300. 
The frequency to is chosen such that hu is constant. The thickness of the PML layer is given by 
Wpmi = 5 except for the case with 300 subdomains which we simulated twice, with u> pm i = 5 and 
Wpml = 6. 

What stands out is that the convergence is very fast, with between 4 and 9 iterations required 
for reduction of the residual by 10 -6 . There is only a mild dependence on the grid size and on 
the number of subdomains. The dependence on w pm i and the somewhat larger number for 300 
subdomains with u> pm i — 5 will be discussed below. 

4-2. Example 2: A random medium in 3-D 

Our second example is random medium in 3-D. Plots of the medium and a solution are given 
in Figure H The size of the example varied between 100 3 and 400 3 (excluding PML layers on the 
sides). In all cases the medium was divided in layers of thickness 10 (excluding again the PML 
layers). Experiments were performed with u> pm i = 4 and 5. The thickness of the subdomain on 
which the computation took place was hence 19 and 21 grid points respectively. The results are 
summarized in Table [3] 

The result are similar to those of the Marmousi examples. The iterative method converged 
rapidly, in 6 to 8 iterations. In these examples the value u> pm i = 4 is sufficient. 
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(a) 

c(x,y) for Marmousi (m/s) 
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(b) 

wavefield for Marmousi (om/(2jt) = 50) 
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Figure 1: Marmousi model and solution with ^- = 50 

4-. 3. Example 3: A constant medium in 3-D and varying w pm \ 

In our third example we explore the dependence of the convergence on u> pm i- The main conclu- 
sion of the previous two examples is that convergence is fast in all cases. Nevertheless, a increase 
in w pm i reduces the number of iterations somewhat in the larger examples. In fact, the iteration 
count stays remarkable in the second column of Table [3] 

In this example the domain is the unit cube, and the velocity c = 1. (A constant medium 
is attractive because it requires less computational resources, due to the fact that for only one 
subdomain the LDL 1 decomposition has to be computed.) The subdomain size varies between 
100 3 and 400 3 (excluding the outer PML layers), while the thickness of the PML layers varies 
between 3 and 6 gridpoints. The frequency us is chosen to correspond to 10 grid points per 
wavelength. The results are given in Table [3] 

While we have limited data, still the following pattern can be observed. For fixed w pm \ the 
number of iterations increases with the grid size. However the number of iterations can be kept 
more or less constant if one can increases if pm i at the same time as the grid size. Here ui pm i goes 
roughly logarithmically with the grid size. 

5. Conclusions 

A new domain decomposition method for the Helmholtz equation was presented. It has re- 
markably fast convergence, even in the case of thin-layered subdomains. We have focussed on the 
use of this method with sparse direct solvers on the subdomains. For this combination there are 
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(a) 




Figure 2: (a) Random medium used in example 2; (b) solution to the Helmholtz equation with a point source. 

still open questions about the best implementation of the method. We like to point the reader 
to [S3] where parallel methods are developed for the moving PML sweeping method, that appear 
applicable in this case also. Our code has not (yet) been optimized this way. Other open questions 
are the applicability for more general domains and more general choices of subdomains, for example 
2-D decompositions of a 2-D domain or 3-D decompositions of a 3-D domain, and the combination 
with other methods for solving the subdomain problems. 

The solutions to the time harmonic Maxwell equations and the time harmonic linear elastic 
wave equation behave in many respects the same as those of the Helmholtz equation. We expect 
that the techniques outlined in this paper are applicable in those cases as well. 
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